Fisk distribution (fisk) — the log-logistic workhorse#

The Fisk distribution (SciPy: scipy.stats.fisk) is widely known in applied statistics as the log-logistic distribution. It is a continuous distribution on the positive real line with Pareto-like tails.

A key representation makes much of its behavior intuitive:

If (X \sim \mathrm{Fisk}(c, \text{scale}=s)) with loc=0, then

[ \log X \sim \mathrm{Logistic}(\mu=\log s,; \sigma = 1/c). ]

So you can think of Fisk as “a logistic distribution on the log-scale”.

What you’ll learn#

  • how to write the PDF/CDF/quantile function (and how SciPy parameterizes them)

  • which moments exist (and why many do not)

  • how to sample with NumPy only via inverse-CDF / logistic tricks

  • how to fit and validate the model with scipy.stats.fisk


import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import optimize, stats

pio.templates.default = "plotly_white"
# CKC convention for Plotly notebooks
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

# Example parameters used throughout
c0 = 5.0
scale0 = 2.0
loc0 = 0.0
Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & Classification#

  • Name: fisk (Fisk / log-logistic distribution; SciPy: scipy.stats.fisk)

  • Type: Continuous

  • Support (standard): (x \in (0, \infty))

  • Parameter space (standard): shape (c>0), scale (s>0)

  • SciPy location/scale: loc \in \mathbb{R}, scale > 0 with [ X = \mathrm{loc} + \mathrm{scale},Z,\quad Z \sim \mathrm{Fisk}(c, 1). ]

We write (2-parameter form):

[ X \sim \mathrm{Fisk}(c, s). ]

The standard form is (\mathrm{Fisk}(c,1)), i.e. stats.fisk(c, loc=0, scale=1).

2) Intuition & Motivation#

2.1 What it models#

The Fisk/log-logistic distribution models strictly positive quantities where:

  • values can vary over orders of magnitude, and

  • the right tail is heavy (polynomial decay), and

  • on the log-scale, the distribution looks approximately logistic.

A nice way to see the logistic link is to rewrite the CDF in terms of (y = \log x) (with loc=0):

[ F(x) = \frac{1}{1 + (s/x)^c} = \sigma\bigl(c(\log x - \log s)\bigr), \qquad \sigma(t) = \frac{1}{1+e^{-t}}. ]

So the “transition” happens around (x \approx s) (because (\log x - \log s \approx 0)).

2.2 Typical real-world use cases#

  • Survival analysis / reliability: time-to-event with a hazard that can rise and fall.

  • Economics: income/wealth-like quantities with heavy tails.

  • Hydrology / environmental extremes: positive heavy-tailed magnitudes.

  • Generative modeling: as a simple positive heavy-tailed component.

2.3 Relations to other distributions#

  • Log-logistic: Fisk is the log-logistic distribution.

  • Logistic on the log-scale: (\log X) is logistic.

  • Burr Type XII: Fisk is a special case of Burr XII (shape parameter (k=1)).

  • Pareto-like tails: (\mathbb{P}(X>x) \sim (s/x)^c) for large (x), so (c) behaves like a tail index.

3) Formal Definition#

We use the common 2-parameter form with shape (c>0) and scale (s>0).

3.1 PDF#

For (x>0):

[ f(x\mid c,s) = \frac{c}{s},\frac{(x/s)^{c-1}}{\bigl(1 + (x/s)^c\bigr)^2}. ]

3.2 CDF#

For (x>0):

[ F(x\mid c,s) = \frac{1}{1 + (s/x)^c} = \frac{(x/s)^c}{1 + (x/s)^c} = \sigma\bigl(c(\log x - \log s)\bigr). ]

The survival function is

[ \bar F(x) = 1 - F(x) = \frac{1}{1 + (x/s)^c}. ]

3.3 Quantile function (inverse CDF)#

For (u\in(0,1)):

[ Q(u) = F^{-1}(u) = s\left(\frac{u}{1-u}\right)^{1/c}. ]

3.4 SciPy parameterization#

SciPy uses:

stats.fisk(c, loc=0, scale=1)

with (z = (x-\mathrm{loc})/\mathrm{scale}) and support (z>0).

def _log1pexp(a: np.ndarray) -> np.ndarray:
    """Stable log(1+exp(a)) for real a (NumPy-only)."""

    a = np.asarray(a, dtype=float)
    out = np.empty_like(a, dtype=float)

    pos = a > 0
    out[pos] = a[pos] + np.log1p(np.exp(-a[pos]))
    out[~pos] = np.log1p(np.exp(a[~pos]))
    return out


def _expit(a: np.ndarray) -> np.ndarray:
    """Stable logistic sigmoid 1/(1+exp(-a)) (NumPy-only)."""

    a = np.asarray(a, dtype=float)
    out = np.empty_like(a, dtype=float)

    pos = a >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-a[pos]))

    ea = np.exp(a[~pos])
    out[~pos] = ea / (1.0 + ea)
    return out


def _logit(p: np.ndarray) -> np.ndarray:
    """Stable log(p/(1-p)) for p in (0,1)."""

    p = np.asarray(p, dtype=float)
    return np.log(p) - np.log1p(-p)


def fisk_logpdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Log-PDF of Fisk(c, loc, scale) in SciPy's parameterization."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.full_like(z, fill_value=-np.inf, dtype=float)
    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    a = c * logz  # log(z^c)

    out[mask] = (
        np.log(c)
        - np.log(scale)
        + (c - 1.0) * logz
        - 2.0 * _log1pexp(a)
    )

    return out


def fisk_pdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """PDF of Fisk(c, loc, scale)."""

    return np.exp(fisk_logpdf(x, c=c, loc=loc, scale=scale))


def fisk_cdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """CDF of Fisk(c, loc, scale)."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.zeros_like(z, dtype=float)
    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    out[mask] = _expit(c * logz)
    return out


def fisk_sf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Survival function 1 - CDF of Fisk(c, loc, scale)."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.ones_like(z, dtype=float)
    out[z <= 0] = 1.0

    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    out[mask] = _expit(-c * logz)
    return out


def fisk_ppf(u: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Quantile function (inverse CDF) of Fisk(c, loc, scale)."""

    u = np.asarray(u, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if np.any((u < 0) | (u > 1)):
        raise ValueError("u must be in [0, 1]")

    out = np.full_like(u, fill_value=np.nan, dtype=float)
    out[u == 0] = loc
    out[u == 1] = np.inf

    mask = (u > 0) & (u < 1)
    if not np.any(mask):
        return out

    out[mask] = loc + scale * np.exp(_logit(u[mask]) / c)
    return out

4) Moments & Properties#

A useful fact is that the right tail is polynomial:

[ \mathbb{P}(X>x) = \bar F(x) = \frac{1}{1 + (x/s)^c} \sim (s/x)^c\quad (x\to\infty). ]

So (c) is a tail index: smaller (c) means heavier tails.

4.1 Raw moments (and existence)#

For (k>0), the raw moment exists iff (c > k), and

[ \mathbb{E}[X^k] = s^k,\Gamma!\left(1+\frac{k}{c}\right)\Gamma!\left(1-\frac{k}{c}\right) = s^k,\frac{\frac{k\pi}{c}}{\sin\left(\frac{k\pi}{c}\right)}. ]

Consequences:

  • Mean exists for (c>1): [\mathbb{E}[X] = s,\frac{\pi/c}{\sin(\pi/c)}.]

  • Variance exists for (c>2): [\mathrm{Var}(X) = s^2\left(\frac{2\pi/c}{\sin(2\pi/c)} - \left(\frac{\pi/c}{\sin(\pi/c)}\right)^2\right).]

  • Skewness exists for (c>3); kurtosis exists for (c>4).

4.2 Location summaries#

  • Median: (\mathrm{median}(X)=s) (and with SciPy loc, the median is loc + scale).

  • Mode (for (c>1)): [\mathrm{mode}(X)=s\left(\frac{c-1}{c+1}\right)^{1/c}.] For (c\le 1), the density is decreasing and the mode is at the lower endpoint.

4.3 MGF / characteristic function#

  • The MGF (M_X(t)=\mathbb{E}[e^{tX}]) does not exist for any (t>0) (polynomial tails cannot beat exponential growth).

  • The characteristic function (\varphi_X(t)=\mathbb{E}[e^{itX}]) exists, but it does not simplify to an elementary closed form; it’s typically computed numerically.

4.4 Entropy#

Using the logistic-on-log-scale representation, the differential entropy is

[ h(X) = 2 + \log\left(\frac{s}{c}\right)\quad \text{(nats)}. ]

(Translation by loc does not change entropy; scaling multiplies it by adding (\log(\mathrm{scale})).)

def fisk_raw_moment(k: float, c: float, scale: float = 1.0) -> float:
    """Raw moment E[X^k] for Fisk(c, scale) (loc=0).

    Returns np.inf if the moment diverges (c <= k).
    """

    k = float(k)
    c = float(c)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")
    if k == 0:
        return 1.0

    if c <= k:
        return float(np.inf)

    a = np.pi * k / c
    return float((scale**k) * (a / np.sin(a)))


def fisk_entropy(c: float, scale: float = 1.0) -> float:
    """Differential entropy in nats for Fisk(c, scale) (loc drops out)."""

    c = float(c)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    return float(2.0 + np.log(scale / c))


def fisk_mode(c: float, loc: float = 0.0, scale: float = 1.0) -> float:
    """Mode of Fisk(c, loc, scale)."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if c <= 1:
        return float(loc)  # density decreases from the boundary

    return float(loc + scale * ((c - 1.0) / (c + 1.0)) ** (1.0 / c))


def fisk_moments(c: float, loc: float = 0.0, scale: float = 1.0) -> dict:
    """Mean/variance/skewness/kurtosis (if they exist) + useful summaries."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    m1 = fisk_raw_moment(1.0, c=c, scale=scale)
    m2 = fisk_raw_moment(2.0, c=c, scale=scale)
    m3 = fisk_raw_moment(3.0, c=c, scale=scale)
    m4 = fisk_raw_moment(4.0, c=c, scale=scale)

    mean = loc + m1 if np.isfinite(m1) else np.inf

    if np.isfinite(m1) and np.isfinite(m2):
        var = m2 - m1**2
    else:
        var = np.inf

    skew = np.nan
    excess_kurt = np.nan

    if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and var > 0:
        mu3 = m3 - 3 * m2 * m1 + 2 * (m1**3)
        skew = mu3 / (var ** 1.5)

    if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and np.isfinite(m4) and var > 0:
        mu4 = m4 - 4 * m3 * m1 + 6 * m2 * (m1**2) - 3 * (m1**4)
        excess_kurt = mu4 / (var**2) - 3.0

    median = loc + scale

    return {
        "mean": float(mean),
        "var": float(var),
        "skew": float(skew) if np.isfinite(skew) else np.nan,
        "excess_kurtosis": float(excess_kurt) if np.isfinite(excess_kurt) else np.nan,
        "median": float(median),
        "mode": fisk_mode(c=c, loc=loc, scale=scale),
        "entropy": fisk_entropy(c=c, scale=scale),
    }


m0 = fisk_moments(c0, loc=loc0, scale=scale0)
m0
{'mean': 2.13791866423119,
 'var': 0.7145293838425228,
 'skew': 2.4852755496867136,
 'excess_kurtosis': 26.556191909249108,
 'median': 2.0,
 'mode': 1.8442158229634555,
 'entropy': 1.083709268125845}

5) Parameter Interpretation#

5.1 Scale s#

  • scale = s sets the median: (\mathrm{median}(X)=s) (when loc=0).

  • On the log-scale, it is a location parameter: (\mu = \log s).

5.2 Shape c#

  • c controls the steepness around the median and the tail heaviness.

  • The tail index is c: (\mathbb{P}(X>x) \sim (s/x)^c).

  • Moments exist only up to order < c.

  • On the log-scale, c is the inverse of the logistic scale: (\sigma = 1/c).

5.3 Hazard shape (survival analysis intuition)#

The hazard rate (h(x)=f(x)/(1-F(x))) (for loc=0) simplifies to

[ h(x) = \frac{c}{s},\frac{(x/s)^{c-1}}{1 + (x/s)^c}. ]

  • For (c \le 1), the hazard is decreasing.

  • For (c > 1), the hazard is typically unimodal (rises, then falls).

def fisk_hazard(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Hazard h(x) = f(x) / (1 - F(x))."""

    f = fisk_pdf(x, c=c, loc=loc, scale=scale)
    s = fisk_sf(x, c=c, loc=loc, scale=scale)

    out = np.full_like(f, fill_value=np.nan, dtype=float)
    mask = s > 0
    out[mask] = f[mask] / s[mask]
    return out


x = np.logspace(-3, 2, 800) * scale0
c_values = [0.7, 1.0, 2.0, 5.0]

fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=["PDF (log x)", "CDF (log x)", "Hazard (log x)"],
)

for c in c_values:
    fig.add_trace(
        go.Scatter(x=x, y=fisk_pdf(x, c=c, scale=scale0), name=f"c={c}", mode="lines"),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=x,
            y=fisk_cdf(x, c=c, scale=scale0),
            name=f"c={c}",
            mode="lines",
            showlegend=False,
        ),
        row=1,
        col=2,
    )
    fig.add_trace(
        go.Scatter(
            x=x,
            y=fisk_hazard(x, c=c, scale=scale0),
            name=f"c={c}",
            mode="lines",
            showlegend=False,
        ),
        row=1,
        col=3,
    )

for j in [1, 2, 3]:
    fig.update_xaxes(type="log", title_text="x", row=1, col=j)

fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_yaxes(title_text="h(x)", row=1, col=3)

fig.update_layout(title="Fisk distribution shape effects (scale fixed)", height=420)
fig.show()

6) Derivations#

6.1 Expectation (raw moments)#

Start from the definition for (k>0):

[ \mathbb{E}[X^k] = \int_0^\infty x^k,\frac{c}{s},\frac{(x/s)^{c-1}}{(1 + (x/s)^c)^2},dx. ]

Use the substitution (t = (x/s)^c), so (x = s,t^{1/c}) and (dx = \frac{s}{c} t^{1/c-1},dt). A pleasant simplification happens:

[ \frac{c}{s}(x/s)^{c-1},dx = \frac{c}{s},t^{(c-1)/c},\frac{s}{c}t^{1/c-1}dt = dt. ]

So

[ \mathbb{E}[X^k] = s^k \int_0^\infty \frac{t^{k/c}}{(1+t)^2},dt. ]

Recognize the Beta-function integral:

[ \int_0^\infty \frac{t^{p-1}}{(1+t)^{m}},dt = B(p, m-p),\quad 0<p<m. ]

Here (m=2) and (p=1+k/c), so the condition becomes (c>k). Then

[ \mathbb{E}[X^k] = s^k,B\left(1+\frac{k}{c},,1-\frac{k}{c}\right) = s^k,\Gamma\left(1+\frac{k}{c}\right)\Gamma\left(1-\frac{k}{c}\right). ]

Using (\Gamma(1+z)\Gamma(1-z)=\frac{\pi z}{\sin(\pi z)}) gives the sine form.

6.2 Variance#

When (c>2),

[ \mathrm{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2, ]

with (\mathbb{E}[X]) requiring (c>1).

6.3 Likelihood (i.i.d. sample)#

For i.i.d. observations (x_1,\dots,x_n > 0) (assuming loc=0), the log-likelihood is

[ \ell(c,s) = \sum_{i=1}^n \log f(x_i\mid c,s) = n\log c - n\log s + (c-1)\sum_i\log\left(\frac{x_i}{s}\right) - 2\sum_i \log\left(1+\left(\frac{x_i}{s}\right)^c\right). ]

There is no closed-form MLE for ((c,s)); we maximize (\ell) numerically. A common trick is to optimize over ((\log c, \log s)) to enforce positivity.

def fisk_loglikelihood(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> float:
    """Total log-likelihood for i.i.d. data under Fisk(c, loc, scale)."""

    return float(np.sum(fisk_logpdf(x, c=c, loc=loc, scale=scale)))


def fisk_mle_loc_fixed(x: np.ndarray, loc: float = 0.0) -> tuple[float, float]:
    """Numerical MLE for (c, scale) with loc fixed.

    Uses a logistic-on-log-scale initialization.
    """

    x = np.asarray(x, dtype=float)
    if np.any(x <= loc):
        raise ValueError("all observations must be > loc")

    x_pos = x - loc

    # Initialization from log-scale logistic approximation
    y = np.log(x_pos)
    mu_init = float(np.median(y))
    q25, q75 = np.quantile(y, [0.25, 0.75])
    iqr = float(q75 - q25)

    # Logistic IQR = 2*s*log(3), with s = 1/c
    c_init = float(max(2.0 * np.log(3.0) / max(iqr, 1e-6), 1e-3))
    scale_init = float(max(np.exp(mu_init), 1e-6))

    def nll(theta: np.ndarray) -> float:
        log_c, log_scale = float(theta[0]), float(theta[1])
        c = float(np.exp(log_c))
        scale = float(np.exp(log_scale))
        return -fisk_loglikelihood(x, c=c, loc=loc, scale=scale)

    res = optimize.minimize(
        nll,
        x0=np.array([np.log(c_init), np.log(scale_init)]),
        method="Nelder-Mead",
        options={"maxiter": 5000},
    )

    log_c_hat, log_scale_hat = res.x
    return float(np.exp(log_c_hat)), float(np.exp(log_scale_hat))


# Quick demo on synthetic data
x_demo = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=2000, random_state=rng)
c_hat, scale_hat = fisk_mle_loc_fixed(x_demo, loc=loc0)

print("true (c, scale) =", (c0, scale0))
print("MLE  (c, scale) =", (c_hat, scale_hat))
true (c, scale) = (5.0, 2.0)
MLE  (c, scale) = (5.081944016481647, 1.9992827793076853)

7) Sampling & Simulation#

Inverse-transform sampling uses the quantile function.

If (U\sim\mathrm{Uniform}(0,1)), then

[ X = Q(U) = \mathrm{loc} + \mathrm{scale},\exp\left(\frac{\mathrm{logit}(U)}{c}\right) = \mathrm{loc} + \mathrm{scale}\left(\frac{U}{1-U}\right)^{1/c}. ]

Algorithm (NumPy-only):

  1. Draw (u_1,\dots,u_n) i.i.d. from Uniform(0,1).

  2. Clip away from exactly 0 or 1 (to avoid ±inf after logit).

  3. Compute z = exp(logit(u) / c) and return loc + scale * z.

Because the distribution is logistic on the log-scale, this is equivalent to sampling

[ Y \sim \mathrm{Logistic}(\mu=\log(\mathrm{scale}), \sigma = 1/c),\qquad X = \mathrm{loc}+e^Y. ]

def fisk_rvs_numpy(
    c: float,
    size: int,
    loc: float = 0.0,
    scale: float = 1.0,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Sample Fisk(c, loc, scale) using NumPy only (inverse CDF)."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if rng is None:
        rng = np.random.default_rng()

    u = rng.random(size)

    # Avoid u=0 or u=1 exactly (would map to +/-inf)
    eps = np.finfo(float).eps
    u = np.clip(u, eps, 1.0 - eps)

    z = np.exp(_logit(u) / c)
    return loc + scale * z


# Monte Carlo validation (choose c>2 so mean/variance exist)
n = 200_000
samples_numpy = fisk_rvs_numpy(c0, size=n, loc=loc0, scale=scale0, rng=rng)

mc_mean = samples_numpy.mean()
mc_var = samples_numpy.var(ddof=0)

print("theory mean", m0["mean"], "MC", mc_mean)
print("theory var ", m0["var"], "MC", mc_var)

# Compare NumPy-only sampler to SciPy sampler (quick 2-sample KS test)
samples_scipy = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=n, random_state=rng)
ks = stats.ks_2samp(samples_numpy[:30_000], samples_scipy[:30_000])
ks
theory mean 2.13791866423119 MC 2.137557641110836
theory var  0.7145293838425228 MC 0.7035045300612408
KstestResult(statistic=0.009700000000000042, pvalue=0.11809773491898767, statistic_location=2.1769456563992517, statistic_sign=-1)

8) Visualization#

We’ll visualize:

  • the PDF with a Monte Carlo histogram overlay

  • the CDF vs an empirical CDF

  • the log-scale view ((\log X) should look logistic)

Because Fisk can be heavy-tailed, a log-x axis is often the clearest.

# PDF + histogram + CDF/ECDF
x_grid = np.logspace(-3, 2, 600) * scale0

# Empirical CDF (subsample for plotting)
sub = np.sort(samples_numpy[:10_000])
ecdf_y = np.arange(1, sub.size + 1) / sub.size

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=["PDF with Monte Carlo histogram", "CDF vs empirical CDF"],
)

# Histogram (density)
fig.add_trace(
    go.Histogram(
        x=samples_numpy,
        histnorm="probability density",
        nbinsx=120,
        name="samples",
        opacity=0.5,
    ),
    row=1,
    col=1,
)

# Theoretical PDF
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fisk_pdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode="lines",
        name="theory pdf",
        line=dict(width=2),
    ),
    row=1,
    col=1,
)

# Theoretical CDF
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fisk_cdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode="lines",
        name="theory cdf",
        line=dict(width=2),
        showlegend=False,
    ),
    row=1,
    col=2,
)

# Empirical CDF
fig.add_trace(
    go.Scatter(
        x=sub,
        y=ecdf_y,
        mode="lines",
        name="empirical cdf",
        showlegend=False,
    ),
    row=1,
    col=2,
)

fig.update_xaxes(type="log", title_text="x", row=1, col=1)
fig.update_xaxes(type="log", title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)

fig.update_layout(height=420, title=f"Fisk(c={c0}, scale={scale0})")
fig.show()


# Log-scale view: log(X) should look logistic
y = np.log(samples_numpy - loc0)
mu_hat = np.median(y)
# logistic scale estimate via IQR: IQR = 2*s*log(3)
q25, q75 = np.quantile(y, [0.25, 0.75])
s_hat = (q75 - q25) / (2 * np.log(3))

# Compare to SciPy logistic on log-scale
logistic_dist = stats.logistic(loc=mu_hat, scale=s_hat)

y_grid = np.linspace(np.quantile(y, 0.01), np.quantile(y, 0.99), 500)

fig2 = go.Figure()
fig2.add_trace(
    go.Histogram(x=y, histnorm="probability density", nbinsx=120, name="log samples", opacity=0.6)
)
fig2.add_trace(go.Scatter(x=y_grid, y=logistic_dist.pdf(y_grid), mode="lines", name="logistic fit"))
fig2.update_layout(
    title="On the log-scale, Fisk becomes logistic",
    xaxis_title="y = log(x)",
    yaxis_title="density",
    height=380,
)
fig2.show()

9) SciPy Integration (scipy.stats.fisk)#

SciPy parameterization:

stats.fisk(c, loc=0, scale=1)
  • c is the shape.

  • loc shifts the support to (loc, ∞).

  • scale rescales (and in the 2-parameter form with loc=0, it equals the median).

Key methods:

  • pdf(x), logpdf(x), cdf(x), ppf(q)

  • rvs(size, random_state=...)

  • fit(data, ...) for MLE

dist = stats.fisk(c0, loc=loc0, scale=scale0)

x_test = np.array([0.5, 1.0, 2.0, 4.0, 8.0])

pdf_scipy = dist.pdf(x_test)
cdf_scipy = dist.cdf(x_test)

pdf_ours = fisk_pdf(x_test, c=c0, loc=loc0, scale=scale0)
cdf_ours = fisk_cdf(x_test, c=c0, loc=loc0, scale=scale0)

print("pdf close?", np.allclose(pdf_scipy, pdf_ours))
print("cdf close?", np.allclose(cdf_scipy, cdf_ours))

# Fitting (MLE) with SciPy; fix loc if you know the support starts at 0.
c_fit, loc_fit, scale_fit = stats.fisk.fit(samples_numpy[:30_000], floc=0)
print("fit (c, loc, scale) =", (c_fit, loc_fit, scale_fit))
pdf close? True
cdf close? True
fit (c, loc, scale) = (4.998041440572122, 0, 2.007897677685366)

10) Statistical Use Cases#

10.1 Hypothesis testing / goodness-of-fit#

Typical workflow:

  1. Fit the distribution parameters (MLE).

  2. Check fit with diagnostics:

    • KS test (quick but less sensitive in the tails)

    • QQ/PP plots (especially important for tail fit)

    • tail-specific checks (log-log plots, exceedance probabilities)

You can also compare candidates (e.g., Fisk vs lognormal vs Weibull) using AIC.

10.2 Bayesian modeling#

There is no conjugate prior for ((c,s)), but Bayesian inference is straightforward with MCMC/VI. A useful modeling trick is to work with (y_i = \log x_i):

[ y_i \sim \mathrm{Logistic}(\mu=\log s,,\sigma=1/c), ]

and add the Jacobian term if you write the model on (x) directly.

Below we’ll do a simple grid posterior approximation with priors on (\log c) and (\log s).

10.3 Generative modeling#

Because sampling is cheap, Fisk is convenient for generating positive heavy-tailed synthetic data. For example, in an accelerated-failure-time style model, you can let the scale depend on covariates:

[ \log X = \beta^\top x + \varepsilon,\quad \varepsilon \sim \mathrm{Logistic}(0, 1/c). ]

# Hypothesis testing + model comparison (AIC) on synthetic data
true_c, true_scale = 2.5, 1.7
x_data = fisk_rvs_numpy(true_c, size=1500, scale=true_scale, rng=rng)

# Fit Fisk
c_hat, loc_hat, scale_hat = stats.fisk.fit(x_data, floc=0)
ll_fisk = float(np.sum(stats.fisk.logpdf(x_data, c_hat, loc=loc_hat, scale=scale_hat)))
aic_fisk = 2 * 2 - 2 * ll_fisk  # 2 params (c, scale) since loc fixed

# Fit lognormal
s_logn, loc_logn, scale_logn = stats.lognorm.fit(x_data, floc=0)
ll_logn = float(np.sum(stats.lognorm.logpdf(x_data, s_logn, loc=loc_logn, scale=scale_logn)))
aic_logn = 2 * 2 - 2 * ll_logn

# Fit Weibull
c_w, loc_w, scale_w = stats.weibull_min.fit(x_data, floc=0)
ll_w = float(np.sum(stats.weibull_min.logpdf(x_data, c_w, loc=loc_w, scale=scale_w)))
aic_w = 2 * 2 - 2 * ll_w

print("Fitted Fisk   (c, scale)", (c_hat, scale_hat), "AIC", aic_fisk)
print("Fitted lognorm(s, scale)", (s_logn, scale_logn), "AIC", aic_logn)
print("Fitted Weibull(c, scale)", (c_w, scale_w), "AIC", aic_w)

# KS test against fitted Fisk CDF (note: p-values are approximate after fitting)
dist_hat = stats.fisk(c_hat, loc=0, scale=scale_hat)
ks = stats.kstest(x_data, dist_hat.cdf)
ks
Fitted Fisk   (c, scale) (2.459344430482952, 1.6658920422291383) AIC 4812.347844437876
Fitted lognorm(s, scale) (0.7263480358234555, 1.661320104043914) AIC 4824.475216282033
Fitted Weibull(c, scale) (1.3345162113961129, 2.3881880476861284) AIC 5114.144353642654
KstestResult(statistic=0.01612416741026612, pvalue=0.8241228923230904, statistic_location=2.6508596780817344, statistic_sign=-1)
# Simple Bayesian grid posterior over (c, scale) with log-normal-ish priors
x_small = x_data[:200]
logx = np.log(x_small)

c_grid = np.linspace(0.8, 6.0, 120)
scale_grid = np.linspace(0.3, 4.0, 140)
log_scale_grid = np.log(scale_grid)

# Priors: log c ~ N(log 2, 0.7^2), log scale ~ N(log 2, 0.7^2)
mu0 = np.log(2.0)
s0 = 0.7

log_prior_c = stats.norm.logpdf(np.log(c_grid), loc=mu0, scale=s0)
log_prior_scale = stats.norm.logpdf(log_scale_grid, loc=mu0, scale=s0)

log_post = np.empty((c_grid.size, scale_grid.size), dtype=float)

for i, c in enumerate(c_grid):
    # Vectorized over scale_grid
    logz = logx[:, None] - log_scale_grid[None, :]
    a = c * logz

    # logpdf for each observation and scale
    logpdf = (
        np.log(c)
        - log_scale_grid[None, :]
        + (c - 1.0) * logz
        - 2.0 * _log1pexp(a)
    )

    loglik = logpdf.sum(axis=0)
    log_post[i, :] = loglik + log_prior_c[i] + log_prior_scale

# Normalize to probabilities
log_post -= np.max(log_post)
post = np.exp(log_post)
post /= post.sum()

# Posterior summaries
c_mean = float((post.sum(axis=1) * c_grid).sum())
scale_mean = float((post.sum(axis=0) * scale_grid).sum())

ij_map = np.unravel_index(np.argmax(post), post.shape)
c_map = float(c_grid[ij_map[0]])
scale_map = float(scale_grid[ij_map[1]])

print("posterior mean (c, scale) =", (c_mean, scale_mean))
print("MAP           (c, scale) =", (c_map, scale_map))

fig = go.Figure(
    data=go.Heatmap(
        x=scale_grid,
        y=c_grid,
        z=post,
        colorscale="Viridis",
        colorbar=dict(title="posterior"),
    )
)
fig.update_layout(
    title="Grid posterior p(c, scale | data) (loc fixed to 0)",
    xaxis_title="scale",
    yaxis_title="c",
    height=420,
)
fig.show()
posterior mean (c, scale) = (2.5669389607324895, 1.697934398969099)
MAP           (c, scale) = (2.547899159663866, 1.6841726618705037)
# Generative modeling sketch: scale depends on covariates (AFT-style)
n_gen = 1500
x_feat = rng.normal(size=n_gen)

beta0, beta1 = 0.2, 0.9
scale_x = np.exp(beta0 + beta1 * x_feat)  # positive, log-linear

u = rng.random(size=n_gen)
u = np.clip(u, np.finfo(float).eps, 1.0 - np.finfo(float).eps)

# Same shape c0, but varying scale per observation
x_gen = scale_x * np.exp(_logit(u) / c0)

fig = go.Figure()
fig.add_trace(
    go.Scatter(x=x_feat, y=x_gen, mode="markers", marker=dict(size=4, opacity=0.5))
)
fig.update_layout(
    title="Synthetic data: Fisk noise with covariate-dependent scale",
    xaxis_title="feature x",
    yaxis_title="generated positive outcome",
    height=380,
)
fig.show()

11) Pitfalls#

  • Invalid parameters: c <= 0 or scale <= 0 is not valid.

  • Moment non-existence:

    • if c <= 1, the mean is infinite/undefined

    • if c <= 2, the variance is infinite

    • more generally, (\mathbb{E}[X^k]) exists only for k < c

  • Parameterization confusion:

    • “Fisk” and “log-logistic” are the same distribution.

    • Many texts use (\alpha) (shape) and (\beta) (scale); SciPy uses c and scale.

    • On the log-scale, the logistic scale is 1/c.

  • Numerical overflow:

    • direct computation of ((x/s)^c) can overflow for large x.

    • prefer logpdf, and compute the CDF via sigmoid(c * log(x/s)).

  • Sampling edge cases:

    • inverse CDF uses logit(u); clip u away from 0 and 1.

  • Fitting with loc free:

    • estimating loc changes the support boundary and can be unstable; fix loc when it is known (e.g., 0 for strictly positive data).

# Numerical stability demo: naive power vs log-space

def fisk_pdf_naive(x: np.ndarray, c: float, scale: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    z = x / scale
    return (c / scale) * (z ** (c - 1.0)) / (1.0 + z**c) ** 2


x_big = np.array([1e2, 1e10, 1e50])
c_test = 5.0

with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
    naive = fisk_pdf_naive(x_big, c=c_test, scale=1.0)

stable = fisk_pdf(x_big, c=c_test, scale=1.0)
log_stable = fisk_logpdf(x_big, c=c_test, scale=1.0)

print("x", x_big)
print("naive  pdf", naive)
print("stable pdf", stable)
print("logpdf    ", log_stable)
x [1.e+02 1.e+10 1.e+50]
naive  pdf [0. 0. 0.]
stable pdf [0. 0. 0.]
logpdf     [ -26.0216 -136.5457 -689.1661]

12) Summary#

  • fisk is the log-logistic distribution: (\log X) is logistic.

  • It models positive heavy-tailed data; the tail index is c.

  • Moments exist only up to order < c (mean requires c>1, variance requires c>2).

  • Sampling is easy via the inverse CDF: x = loc + scale * exp(logit(u)/c).

  • SciPy provides a robust implementation: scipy.stats.fisk (pdf/cdf/rvs/fit).

References

  • SciPy: scipy.stats.fisk

  • Burr Type XII distribution (Fisk as a special case)

  • Logistic distribution entropy and quantiles (for the log-scale view)